Systematic simulations of modified gravity: symmetron and dilaton models 
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We study the linear and nonlinear structure formation in the dilaton and symmetron models of modified 
gravity using a generic parameterisation which describes a large class of scenarios using only a few parameters, 
such as the coupling between the scalar field and the matter, and the range of the scalar force on very large scales. 
For this we have modified the iV-body simulation code ECOSMOG, which is a variant of RAMSES working in 
modified gravity scenarios, to perform a set of 110 simulations for different models and parameter values, 
including the default ACDM. These simulations enable us to explore a large portion of the parameter space. We 
have studied the effects of modified gravity on the matter power spectrum and mass function, and found a rich 
and interesting phenomenology where the difference with the ACDM template cannot be reproduced by a linear 
analysis even on scales as large as k ~ 0.05 /iMpc - . Our results show the full effect of screening on nonlinear 
structure formation and the associated deviation from ACDM. We also investigate how differences in the force 
mediated by the scalar field in modified gravity models lead to qualitatively different features for the nonlinear 
power spectrum and the halo mass function, and how varying the individual model parameters changes these 
observables. The differences are particularly large in the nonlinear power spectra whose shapes for f(R), dilaton 
and symmetron models vary greatly, and where the characteristic bump around 1 /iMpc^ 1 of f(R) models is 
preserved for symmetrons, whereas an increase on much smaller scales is particular to symmetrons. No bump is 
present for dilatons where a flattening of the power spectrum takes place on small scales. These deviations from 
ACDM and the differences between modified gravity models, such as dilatons and symmetrons, could be tested 
with future surveys. 



I. INTRODUCTION 

The apparent acceleration of the Universe could be due to at 
least four different reasons: a cosmological constant, dark en- 
ergy 01, modified gravity Q or large spatial inhomogeneities 
y|]. The last of these violates the Copernican principle and re- 
quires a theory for the initial conditions of the Universe while 
the first three invoke a change of the dynamics of the Universe 
itself. 

The cosmological constant solution is rather peculiar as no 
real dynamics is attached to it until the vacuum energy starts 
dominating the energy content of the Universe. This seems to 
have happened in the quite recent past, a fact which is prob- 
lematic and related to the astoundingly small value of the crit- 
ical density of the Universe compared to particle physics ex- 
pectations, which scale as the fourth power of the mass of any 
heavy particle present in the early Universe. 

To alleviate this problem, two other possibilities are com- 
monly invoked. The first one is dark energy [Q]], in which the 
dynamics of a field (e.g., a scalar field in the simplest case) de- 
termines the fate of the Universe. So far no real solution to the 
cosmological constant problem has been found within this set- 
ting although phenomenological works abound. Setting aside 



the problem of the actual value of the dark energy density 
now, these models suffer from another serious problem: dark 
energy evolves on cosmological time scales only when the 
scalar field leads to a long range interaction. Of course, one 
can decree that dark energy does not couple to baryons as in 
coupled quintessence modelfl and therefore alleviate gravita- 
tional problems linked to the existence of a scalar fifth force. 
If this is not the case, then a solution which has been put for- 
ward in the last decade is screened modified gravity mediated 
by a scalar field. 

Many models of screened modified gravity have been con- 
structed so far, which fall within two broad categories. Fol- 
lowing the initial works on massive gravity, models involv- 
ing nonlinear kinetic terms, such as the Galileon |4j-|6|], make 
use of the Vainshtein mechanism whereby large nonlin- 
earities in the vicinity of dense objects effectively reduce the 
scalar coupling to matter to be below the experimental bounds. 
Another class of models originating from the chameleon the- 
ory iH lit] use a screening of the fifth force in dense environ- 
ments due to the nonlinearities of either the scalar potential 
or its coupling to matter (or both). Chameleon models such 
as f(R) gravity lflol - H2ll are such that the mass of the scalar 
field becomes large in dense bodies, effectively suppressing 
the magnitude of the scalar force; other models such as the 
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1 We regard the coupled quintessence model as an example of dark energy 
rather than modified gravity, for which we require a universal coupling to 
all matter species. 
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dilatons lfl3ll and symmetrons lfl4j[l5ll are such that the effec- 
tive coupling to matter becomes vanishingly small in dense 
environments. All cases in the second class of screened mod- 
ified gravity can be described by the same formalism which 
has been recently unified lfl^[l7ll . In this paper, we will con- 
centrate on the second class. 

It has been shown in lfl7ll that the background cosmology 
of these models is extremely constrained. Indeed, the fact that 
particle masses (in the Einstein frame) and the gravitational 
constant (in the Jordan frame) cannot vary substantially be- 
tween the era of the Big Bang Nucleosynthesis (BBN) and 
now implies that the scalar field must stay very close to the 
minimum of the effective potential since before BBN. This is 
guaranteed when the mass of the scalar field on the cosmolog- 
ical background is much heavier than the Hubble expansion 
rate, securing the stability of the minimum to 'kicks' occur- 
ring when particles such as the electrons decouple fl8ll . A 
consequence of this is that the effective equation of state of 
the scalar field in the late-time Universe becomes extremely 
close to —1, hardly distinguishable from the pure A-cold dark 
matter (ACDM) scenario. In practice, models of f(R) gravity, 
chameleon, dilaton and symmetron types usually behave like 
ACDM in the background cosmology since before BBN. 

Fortunately, this does not imply that their cosmology is to- 
tally degenerate with that of the ACDM model: the effects 
of modified gravity appear in the structure formation. Indeed, 
within the Compton wavelength of the scalar fielcfl gravity is 
modified and the growth rate of structures is altered [fnl H . 
At the linear level, this results in a modification of the growth 
equation which depends on the scalar field mass m(a) and the 
coupling to matter (3(a) expressed as functions of the scale 
factor. It turns out that all screened modified gravity models 
with no higher derivative terms in their Lagrangian, including 
their field-dependent potential V(ip) and the coupling to mat- 
ter (3(ip), can be fully reconstructed from the sole knowledge 
of the functions to (a) and /3(a). This allows one to engineer 
models directly from their linear perturbation properties, i.e., 
given m{a) and /3(a) one can build a fully consistent model 
of modified gravity defined by P(ip) and V(ip) llT6[[l7ll . which 
implies that one could study the nonlinear evolution of cos- 
mic structures in the late Universe simply from the knowl- 
edge of m(a) and /3(a). This provides a systematic approach 
to screened modified gravity which can be applied to gener- 
alised chameleon, dilaton and symmetron models. For other 
schemes to parameterise modified gravity see lfl9l - l24ll . 

Studying the nonlinear regime of structure formation is of 
particular importance for screened modified gravity models , 
as local gravity tests often imply that deviations from gen- 
eral relativity are strongest on megaparsec (Mpc) scales IU7I1 . 
where nonlinearities cannot be neglected. Two competing ef- 
fects influence the dynamics of modified gravity here. On the 
one hand, the gravitational interaction is enhanced by the pres- 
ence of a long-range fifth force which implies an increase of 



2 The Compton wavelength of a scalar field is defined as A = m of ^ , and 
m e ff is the effective mass of the scalar field (see below). 



the growth of structure. On the other hand, where local mat- 
ter densities are high enough, screening effects develop and 
structure formation converges to its GR behaviour. These two 
competing effects have been confirmed in already-available 
A-body simulations of f(R) gravity lf25l-[34ll . chameleon ll35l - 
O, dilaton Q and symmetron <4lU42ll models. 

In this work, we apply the (m(a), /3(a)) parameterisation 
to generalise dilaton and symmetron models and study their 
large-scale structure formation. We use modified versions of 
the ECOSMOG code 0] to run A-body simulations in these 
models. This code is based on the publicly-available adaptive 
mesh refinement (AMR) code RAMSES Q, which is effi- 
ciently parallelised and suitable to run simulations systemati- 
cally. The AMR nature of the code means that a higher reso- 
lution can be achieved, without sacrificing the overall perfor- 
mance of the code, in dense regions where the field equations 
are most nonlinear, ensuring the accuracy of the fifth force cal- 
culation there. As a result, our simulations are able to probe 
the structure formation in these modified gravity models down 
to scales well below the typical dark matter halo sizes. 

The results of our simulations indicate that large deviations 
from ACDM in the power spectrum can be found on scales of 
order 1 Mpc for both symmetron and dilaton models for val- 
ues of the parameters which comply with the local constraints 
(the gravitational tests in the Solar system and a mild sup- 
pression of the fifth force on galactic scales typically impose 
that the range of the fifth force should be less than a few Mpc 
in the cosmological background). Large differences are also 
present in the number density of intermediate-sized dark mat- 
ter halos with masses of order 10 13 — lO 14 /!" 1 AI@ (represent- 
ing objects from groups of galaxies to small galaxy clusters). 
For models with a fifth force whose range in the cosmological 
background is of order Mpc and a coupling strength to matter 
of order unity, the deviation from ACDM can reach ~ 40% 
in the symmetron case and ~ 30% in the dilatonic one. Such 
large differences are testable using future galaxy surveys. 

Moreover, symmetron and dilaton models are distinguish- 
able thanks to the very different time dependence of their cou- 
plings to matter. For symmetrons, the coupling has a slow de- 
pendence on the scale factor a in the recent past of the Uni- 
verse and vanishes before a transition redshift (its defini- 
tion will be given later). Dilaton models have a much sharper 
dependence on the scale factor and generically decrease ex- 
ponentially fast going back in time. As will be discussed in 
detail in § III BL the time dependence of the coupling strength 
can be roughly translated into a density dependence, and the 
steep density dependence in the recent past of the Universe (or 
equivalently in regions of low matter density) for dilaton mod- 
els suggests that the dilaton screening is more efficient. These 
properties make the matter power spectra and halo mass func- 
tions behave qualitatively differently in these models. We will 
give a more detailed summary of the results in the concluding 
section. 

The layout of this paper is as follows: in § [II] we review 
scalar-tensor theories and show how such theories of modified 
gravity can be analysed using a simple parametrisation which 
encapsulates all the dynamics; in § [III] we briefly describe the 
generalised symmetron (§ MI Al l and dilaton (§ IIIIBI ) mod- 
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els and the possible effects of varying each model parameter; 
the equations that will be used in the iV-body simulations are 
summarised in § II V 1 while the details are given in § |1V B| we 
next carry out tests of our codes in § [V] and the cosmologi- 
cal simulations of this work are then discussed in § |VI]for the 
symmetron (§ I VI Al l and dilaton (§ I VI Hb cases respectively; 
finally we summarise and conclude in § I VIII 

In the paper we use the units H = c = 1 except where c ap- 
pears explicitly. Overbar (subscript o) denotes the background 
(present-day) value of a quantity and subscript v means 
djdip. k = 87tGn = , where Mp\ is the reduced Planck 
mass and Gat is Newton's constant, are used interchangeably. 



has a minimum <p m i n (Pm)- The mass of the scalar field at the 
minimum, 



d 2 V cS 
dip 2 



(7) 



must be positive. In a cosmological setting we will also im- 
pose that m 2 ^> H 2 with H being the Hubble expansion rate. 
This guarantees the stability of the minimum to perturbations. 
When matter is described by a pressure-less fluid with 



= Pm u»u" 



(8) 



where u M = dx^/dr is the 4-velocity field of the fluid and r 
is the proper time, the matter density p m is conserved 



n. MODIFYING GRAVITY WITH A SCALAR FIELD 

A. Screened modified gravity 

The action governing the dynamics of a scalar field ip in a 
scalar-tensor theory is of the general form 



S 



d^xJ^ 



9 



d A x 



(1) 



where g is the determinant of the metric g M „, R is the Ricci 

scalar and ip m are various matter fields labelled by i. A key 
ingredient of the model is the conformal coupling of <p with 
matter particles. More precisely, the excitations of each mat- 
ter field ipm couple to a metric g^ u which is related to the 
Einstein-frame metric g^ by the conformal rescaling 

9u.v = A 2 {tp)g liV . 



(2) 



The metric g M „ is the Jordan-frame metric. The fact that the 
scalar field couples to matter implies that the scalar field equa- 
tion becomes density-dependent. More specifically, the scalar 
field equation of motion (EOM) is modified due to the cou- 
pling of the scalar field ip to matter: 



Dip = -j3T- 



av 

dip ' 



(3) 



where T is the trace of the energy momentum tensor T^, 
□ = V M V M and the coupling of ip to matter is defined by 



/%) = M P1 



din ,4 

dip 



(4) 



This is equivalent to the usual scalar field EOM with the ef- 
fective potential 



VM = V( V )-[A(<p)-l] T. 



(5) 



We will always require that the effective potential possesses a 
unique density-dependent minimum in the presence of pres- 
sureless matter for which T = —p m , i.e., that the potential 



V eS {<p) = V(ip) + [A(ip) - IK 



(6) 



Pm + Opm = 0, 



(9) 



where 9 = V M u M = ZH is the expansion scalar and the tra- 
jectories are determined by the modified geodesies 



P-^-u 



-0 



V^ip 



Mpi Mpi 
In the weak-field limit with a line element 



-(1 + 2(j))At 2 + (1 - 2cj))dx l dx t 



(10) 



(11) 



and in the non-relativistic case, this reduces to the modified 
geodesic equation for matter particles 



-V 1 [</> + hi A{(p)]. 



(12) 



This can be interpreted as the motion of a particle in the effec- 
tive gravitational potential defined as 



* = </> + lnA(y>), 



(13) 



and is a manifestation of the dynamics of modified gravity. 
One may also call the deviation from the Newtonian gravity a 
fifth force. In this paper we will use these terminologies inter- 
changeably. 

When a particle of mass M in a homogeneous background 
matter density is the source of gravity, the scalar field satisfies 



m 2 ) ip 



(14) 



in which 5^ (r) is the 3-dimensional Dirac <5-function and m 
the scalar field mass in the background. This implies that 



* = -(! + 2/3 2 e- 



G N M 



(15) 



When j3 ~ 0(1) and mr < 1, this implies a substantial devi- 
ation from Newton's law. For bodies much bigger than a point 
particle, nonlinear effects imply that the effective coupling felt 
by a test mass near the source can be much smaller than 1 or 
the scalar field mass becomes much larger than the inverse of 
the typical size of the source (m^ 1 <C r). The dilaton and 
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symmetron models satisfy the first criterion which guarantees 
that solar system and laboratory tests of gravity are evaded. 

In addition to the se//-screening described above, the mod- 
ification of gravity depends on the environment of the bodies 
as well. For example, in a high-density background, the scalar 
field mass m in Eq. ( fl~4T > can be very large, which suppresses 
the deviation from Newtonian gravity according to Eq. ([T3V 

This environmental dependence is at the heart of the screen- 
ing mechanisms in chameleon, dilaton and symmetron cases. 
Indeed, as shown in lfl7tl . the screening is effective when the 
Newtonian potential $ jy generated at the surface of a dense 
body satisfies 



\p 2 



<p c \ < 2/3 00 Mpi$ 



N, 



(16) 



where ip c ,oo are respectively the minimum of the effective po- 
tential inside and far away from the dense body; $n is the 
Newton potential at the surface of the body and (3^ = (3((poo) 
is the coupling to matter outside. Note that the self and envi- 
ronmental screenings are encoded in <5>at and ipoo , (3^ respec- 
tively 

In cosmological simulations, tp^ = p is the background 
value of cp, while (p c is the value inside clustered structures, 
which can be very small. In general, ip could change by several 
orders of magnitude from low-density to high-density regions, 
and this is why the accurate calculation of ip is a challenging 
task. The equations of motion which govern the dynamics of 
the modified gravity models which we consider here are 



V 2 (/> ps 4ttG (p m - p m ) , 
c 2 V 2 <p ps V v {<p) - V v (Cp) + A v (tp) Prl 



. dr 



(17) 

A v (cp)p m ,(l8) 
(19) 



— = -V0 - c 2 p{<p)Vp - P{<p)<P^, 

where in Eq. (fT7T - IT8b we have worked in the quasi-static limit 
so that terms involving time derivatives have been dropped; 
this is a good approximation throughout the course of cosmic 
evolution as the time derivatives are generally much smaller 
than the spatial onefl The first of these equations is the Pois- 
son equation while the last one is the modified Newtonian dy- 
namics due to the presence of the scalar field <p, c.f. Eq. ( fTOb . 
We have reinstated the factors of c because in code units (see 
below) c is no longer unity. 



B. Tomography 

We shall always consider the cosmological evolution of the 
scalar field ip in modified gravity models with a minimum of 
V e ff (ip) at which the scalar field mass m satisfies m 2 3> H 2 . 



3 This has been shown explicitly in, e.g., 1251 . which compares the two di- 
rectly. A more rigorous proof of the validity of the quasi-static approxima- 
tion would be by solving the full time-dependent scalar field EOM, which 
is beyond the scope of the current work. However we find that, in the lin- 
ear perturbation calculations of [17], one gets indistinguishable results by 
solving the full (linearised) EOM and using the quasi-static approximation, 
showing that the latter is actually quite reasonable. 



The time evolution of the scalar field is tightly constrained 
by BBN physics due to its coupling to matter particles. The 
fact that the scalar field evolves along the minimum of V e g(ip) 
implies that the masses of fundamental particles 



m.0 = A(ip)m ha 



(20) 



in which mbare is the bare mass appearing in the matter La- 
grangian, evolve too. In practice, tight constraints on the time 
variation of masses since the time of BBN 



Am,/, Aip 



Mi 



(21) 



pi 



where Atp is the total variation of the field since BBN, impose 
that Am^ /m^f, must be less than ~ 10%. At a redshift of or- 
der z e ps 10 9 , electrons decouple and give a 'kick' lfl8ll to the 
scalar field which would lead to a large violation of the BBN 
bound. To avoid this, the field must be close to the minimum 
of V e s(ip) before z e and simply follow the time evolution of 
the minimum. Moreover, the total excursion of the scalar field 
following the minimum must be small enough. In practice, we 
will always assume that | ip /Mpi | <C 1 along the minimum tra- 
jectory, implying that the BBN bound for the time dependent 
minimum is always satisfied. The models are then valid pro- 
vided the electron 'kick' does not perturb the minimum too 
much. The minimum of the effective potential acts as a slowly 
varying cosmological constant. Indeed, when m 2 ^> H 2 the 
minimum is stable for all the models we will consider. In this 
case, the dynamics are completely determined by the mini- 
mum equation 



dV 

dtp 



-PA 



Mpi 



(22) 



In fact, the knowledge of the time evolution of the mass m 
and the coupling (3 is enough to determine the time evolution 
of the field. Using the minimum equation, we can deduce that 
the field evolves according to 



dip 

~dl 



Mi 



(23) 



pi 



This is the time evolution of the scalar field at the background 
level since the instant when the field starts being at the min- 
imum of the effective potential. The knowledge of the time 
evolution of the mass m and the coupling f3 is enough to deter- 
mine the bare potential V(ip) and the coupling function A(ip) 
completely. To see this, integrating Eq. (|23l once, we find 



cp(a) 



M- 



pi 



am? (a) 



p m (a)da + tp c 



(24) 



where ip c is the initial value of the scalar field at a ull < cibbn 
and we have taken A(tp) ps 1 given that the temporal variation 
of fermion masses must be very weak. If the coupling strength 
j3 is expressed in terms of the field ip and not the scale factor 
a, this is also equivalent to 



1 



dip 3 

I3((p) Mpi 7 a . n . am 2 {a) 



-p m (a)da. 



(25) 
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Similarly the minimum equation implies that the potential can 
be reconstructed as a function of time 



V = V 



M 2 



/? 2 (g) 

am 2 (a) 



p 2 m (a)da, 



(26) 



where Vq is the value of the potential at a = aim- This de- 
fines the bare scalar field potential V(ip) parametrically when 
/3(a) and m(a) are given. Hence we have found that the full 
nonlinear dynamics of the theory can be recovered from the 
knowledge of the time evolutions of the mass and the coupling 
to matter since before BBN. 

The reconstruction mapping gives a one-to-one correspon- 
dence between the scale factor a and the value of the field 
ip(a) in the cosmic background. As the scale factor is in a one- 
to-one correspondence with the matter energy density p m (a), 
we have obtained a mapping p m —> <p(p m ) defined using the 
time evolution of m(a) and /3(a) only. Given these evolutions, 
one can reconstruct] the dynamics of the scalar field for den- 
sities ranging from cosmological to solar system values using 
Eq. d24l i and Eq. (F26l i. By the same token, V(ip) can be recon- 
structed for all values of <p (and p m ) of interest, from the solar 
system and Earth to the cosmological background today. 

In particular, we can now state the screening condition of 
modified gravity models [c.f. Eq. ( TToT ll as 



am 2 (a) 



(27) 



with constant matter densities Pm.out = Pm(» = ftin.out) in- 
side and outside the dense body respectively, and where we 
have defined /3 out = j3(a — a out ). Note that the gravitational 
properties of the screened modified gravity models can be cap- 
tured by the cosmological evolutions of the scalar field mass 
and coupling function only. 

The loosest screening condition follows from the fact the 
Milky Way should be screened as otherwise large deviations 
from Newtonian gravity would have been detected in the solar 
system. For the Milky Way, the density is around six orders of 
magnitude larger than the cosmological background implying 
that ain ~ 10 -2 ; its Newtonian potential is $g ~ 10~ 6 . Tak- 
ing the outside environment to be close to the cosmological 
background we have a out ~ 1. Writing 



m(a) = m /(a), /3(a) = fi g(a), 



(28) 



where / and g are smooth functions of a with slow variations 
we find 



3n m0 H 2 f 1 g(a) 



a,„ « 4 / 2 («) 



-da < M^ G , 



(29) 



in which f2 m is the fractional matter density. Defining I = 

f 1 fro/ \ da, we find that 



3f2 m n7 



— a- > - 

h 2 ~ a> 



(30) 



4 This is done by assuming that the scalar field always minimises its effective 
potential V c ff , and thus the results below are more of qualitative estimates 
than quantitatively accurate predictions. 



Typically this implies that mo/Ho > 10 3 . Hence we find that 
screened models of modified gravity can only act on scales 
below the order of a few Mpc. In fact we will make use of the 
ratio 



i 



to ' 



which is related to the range of the fifth force as 

A = 2998£ h^Mpc. 



(31) 



(32) 



These scales, in the Mpc range, are beyond the linear pertur- 
bation regime and can only be accurately analysed using nu- 
merical simulations. This is the aim of the present article. In 
the next subsection, we will describe the models we will study 
in detail numerically. 



C. The dilaton and symmetron models 

1. Dilatons 

The environment-dependent dilaton model was originally 
described in |[13m . The essential features of the dilaton model 
include a runaway potential and a coupling function A((p) 
which has a minimum. The potential is derived in the strong 
coupling limit of string theory and the form of the coupling 
function ensures the field does not runaway to infinity, which 
would imply decompactification. In lfl3ll the coupling function 
and bare potential of the scalar field were specified as follows: 



M<P) = l+l^-{<P-<P.f 
z ju p1 

V(ip) =V e-™ /Mpi . 



(33) 
(34) 



Here A2 ^> 1, 7 > are dimensionless model parameters, Vq 
is a model parameter with mass dimension 4 and <p* an arbi- 
trary constant. The screening mechanism of the dilaton model 
is shown in Fig.[TJ Again, denoting the value of ip which min- 
imises V e s(f) by ip m in, when matter density is high (p m i n is 
very close to ip* so that f3((p m i n ) w /3(<£*) = and the fifth 
force essentially vanishes, while when matter density is low 
tp m in can evolve away from ip* so that f3(ip m i n ) ^ P(<p*) = 0' 
giving rise to a non-negligible fifth force. 

To study the cosmology of the dilaton model we need only 
consider the dynamics in the vicinity of the field 93*, where 



A 2 
Mpi 



from which we deduce that 



In 



da 



and therefore 

\P(<p)\ = |/% c )|exp 



9A 2 f2 TO o#o / — 1 — tt\ 
Ja lni a 4 m 2 (a) 



(35) 



(36) 



(37) 
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FIG. 1. An illustration of how the dilaton mechanism works. The dashed, dotted and solid curves are respectively the bare potential V((p) of the 
dilaton field, the coupling function and the total effective potential V B a (<p). Left Panel: in high matter-density regions the minimum of V e B (ip) 
is where the coupling strength vanishes and so the fifth force is suppressed. Right Panel: in low matter-density regions the coupling strength 
does not vanish at the minima of V c g (ip), where the dilaton field resides, and so a nonzero fifth force takes effect in structure formation. 




FIG. 2. An illustration of how the symmetron mechanism works. The dashed, dotted and solid curves are respectively the bare potential V(p>) 
of the symmetron field, the coupling function and the total effective potential V e ff (</>)• Left Panel: in high matter-density regions the minimum 
of V c g(tp) is where the coupling strength vanishes and so the fifth force is suppressed. Right Panel: in low matter-density regions the coupling 
strength does not vanish at the minima of V e ff (</?)> where the symmetron field resides, so a nonzero fifth force takes effect in the structure 
formation. 



This is the relation between the coupling at the initial time and 
other cosmological times. 



The initial coupling (taken at a- m \ < obbn) is the same as 
in dense matter on Earth and is related to the cosmological 



value of (3 today, (3{<po), by 



|/% )| = |/%c)|exp 



do 



am 2 (a) 



(38) 



It is possible to have a very small coupling in dense matter 
(|/3((/? c )| <C 1) for any value of the coupling on cosmological 
scales (|/3(iy9o)|) provided that Ai > and that the time varia- 
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tion of m(a) is slow and does not compensate the 1/a 4 diver- 
gence in the integrand. In this situation, the coupling function 
(3 converges exponentially towards zero: this is the Damour- 
Polyakov mechanism 14011 . The fact that A-2 > guarantees 
that the minimum of the coupling function A(<p) is stable and 
becomes the minimum of the effective potential which attracts 
the scalar field at late times. If A2 < 0, the effect of the 
coupling is destabilising and implies that p> diverges exponen- 
tially fast away from tp*. 

Alternatively, a smooth variation of the coupling function to 
matter in the cosmological background and therefore interest- 
ing consequences for the large-scale structure can be achieved 
when the evolution of the mass of the scalar field compensates 
the 1/a 4 factor in the radiation era and evolves in the matter 
era. This is obtained for models with 

m 2 (a) = 3A 2 H 2 (a)M$ 1 . (39) 

Indeed, H(a) ~ aT 2 in the radiation era, which implies that 
the time variation of (3(ip) between BBN and matter-radiation 
equality is 



P(tp) = (3{ip c ) exp 



(40) 



in which Q, r is the fractional density for radiation, and in the 
matter-dominated era 



'<(i 



(41) 



in which a subscript cq denotes the value of a quantity at the 
matter-radiation equality. This is the behaviour of the dilaton 
models already analysed in ll39Tl . 



2. Symmetron 

The symmetron model was originally described in Ill4l 
for which the coupling function and bare potential of the 
scalar field take the following forms respectively: 



V(<p) =V - i/iV + Ja/. 



(42) 



(43) 



Here M < 10 _3 Mpi is a mass scale and /x ~ Ho, A <C 1 
are model parameters. The screening mechanism of the sym- 
metron model is shown in Fig. [2] When the matter density 
is high cpmin coincides with the minimum of A(ip) such that 
P(<Pmin) = and the fifth force vanishes, whilst when mat- 
ter density is low /3(p m i n ) 7^ 0, resulting in a cosmologically 
interesting fifth force. 

A fundamental property of the symmetron models is that 
the coupling to matter vanishes identically in dense regions or 
at redshifts z > z*, and an order-unity coupling is obtained 
after a transition at a redshift z* and in the low matter-density 
regions. In the original symmetron model, this is given by 



P(a) = M 1 



(44) 



for z < z* and (3 = for z > z*. Similarly, 

m(a) = to* 



(45) 



Notice that for symmetron models a subscript * denotes the 
value at far future (a — > 00), and a subscript » means the value 
at the symmetry breaking, i.e., when /3(a) becomes nonzero in 
the cosmological background. 

Using the reconstruction mapping, it is straightforward to 
find that 



rt«) = vWi-(£) 



for z < z* and ip = before. Here we have defined 



and 



1* = V2> 



V, P* 



'Pl 



PmO a * 



(46) 



(47) 



(48) 



The potential for z < z» as a function of a can then be recon- 
structed, using the technique introduced above, as 



V(a) =V + 



Pipi 



2to2M|, 



/a* 
V a 



- 1 



(49) 



The potential as a function of tp can then be found to take the 
form of Eq. (l43l , with fi given in Eq. d48b and 

, P 2 

A = ^. (50) 
Meanwhile, (3 as a function of <p is reconstructed as 

Pip) - —<P- (5D 
<P* 

It could be checked that this agrees with Eq. (l42l . by taking 
(3 = din A/ dip w dA/dp, where the w symbol comes from 
the fact that Awl. 



III. GENERALISED SYMMETRON AND DILATON 
MODELS 

In this section we discuss the generalisations of the dilaton 
and symmetron models, and the effects of varying the model 
parameters. 



A. Generalised symmetron model 

1. Model parameterisation 

The original symmetron model discussed in the previous 
section only includes one specific potential. As a straightfor- 
ward generalisation of this idea, let us consider the following 
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m(a) and (3(a): 



m(a) 



/3(«) = A 



1 - 




(t)1 







a 



(52) 
(53) 



where m, h are two new parameters and not necessarily equal 
to each other, and (m+ , 8+) are the mass and coupling in vac- 
uum as above. As in IU6tl . if the scalar field always follows^ 
^min, one can obtain the following solution for ip(a): 



a / 
3 



n — 2m-f 1 



(54) 



where we have defined tp+ = fl _ 2 ^ t+1 fl m 8+t; 2 a~ 3 and from 
here we will neglect the subscript o in f2 m o. Note that Eq. ( T54l 
is only valid if n — 2m +1^0; the case of h — 2m = — 1 
corresponds to a potential that is not bounded below and is 
therefore not a viable physical model. Again, Eq. (154-b is for 
a > a* and for a < a* we have <p(a) = 0. 

To study the nonlinear evolution of cp, we need to know 
V v (tp) as it appears in the A-body equations Eq. ( fT8l . Noting 
that tp increases monotonically with a, we find 



V, 



_ d[V(a)] da 
da dip 

= — (h — 2m + l)m+ (p* 
= — (h — 2m + l)m+ip* 



1- - 



1-1^- 

p* 



Defining the parameters 

2n - 2m + 2 

M 



N 



2n - 2m + 1 



n — 2m +1 n — 2m + 1 

we find that the potential can be written quite simply as 



V{<p) 



£ 2 (M - N) 



1 

'N 



/ \ N 

(i) < 













In a similar manner, for a > a* we get 

/%) = B(a(<p)) = 8* 



N-l 



(55) 



(56) 



•(57) 



(58) 



It is evident that when N = 2 and M = 4 we recover the orig- 
inal symmetron model. In what follows we will only consider 
M, N to be even and positive integers with M > N to avoid 
having a potential that is unbounded from below. 



2. Effects of varying model parameters 

Let us analyse the effects of varying the five model param- 
eters a*, B+, N, M and £ on structure formation. 

As discussed in lfl7tl . the modifications of the structure for- 
mation at the linear perturbation level is completely deter- 
mined by the two temporal functions m(a) and 8(a), from 
which we can see that: 

1. The strength of the fifth force vanishes for a < a* and 
approaches 28 2 times that of the Newtonian gravity for 
a ^> a*. Decreasing a* increases the time during which 
the fifth force is active thus enhances the matter cluster- 
ing today. 

2. Increasing 8* makes 8 larger at all times, which makes 
the fifth force stronger and leads to more clustering. 



' See 1 50] for a more detailed discussion on the time-evolution of tp. 



3. According to Eq. ( T58l . increasing N makes 8 smaller 
because \ip\ < \ip+\ in general. This can weaken the 
effect of the fifth force. It is because of this reason 
that the symmetron screening is more efficient than the 
chameleon screening with a constant 8 Il7ll . 

4. By increasing M the scalar field will make the transi- 
tion from tp = to >p = ip+ much quicker, because then 
<p M is smaller for small tp and so (1) the symmetry in 
V e s(<p) is easier to be broken and (2) V c a(tp) becomes 
steeper from tp — to p> = p>± . This leads to a stronger 
(and earlier kick-in of the) fifth force and thus matter 
becomes more clustered. 

5. An increase in £ is equivalent to an increase in the range 
A* of the fifth force since A* = 2998£ Mpc/h in vac- 
uum. This extends the modifications of gravity to larger 
cosmological scales and decreases the exponential fac- 
tor e~ mr of suppression of the fifth force. 

These properties will be investigated in depth using A-body 
simulations below. 



B. Generalised dilaton model 

1. Model parameterisation 

The environment-dependent dilaton model has already been 
presented in the previous section. For the model in lH3ll it can 
be shown that 

m(a) = TOoa _ 2 ; (59) 

8(a) = 8 Q a 9n - A ^\ (60) 

where both m(a) and 8(a) are power law functions of a. If 

m(a) = m^a~ r , (61) 

with r 7^ 3/2, then 8 is no longer a power law function of a, 
as we will see below. 

As a straightforward generalisation of the dilaton idea, let 
us consider a quadratic coupling function A(tp) which has a 
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minimum at tp*. Near tp* we have /3(tp) « A2(p — ip*)/Mp\. 
Assuming that the dilaton field always follows the minimum 
of V e s(p), p m - m , one can solve for /3(a) from an integral Il6ll : 



(a < 1) = /3 exp 
= /3 exp 



2r-4 



da 



2r-3 



(a 



2r-3 



1) 



(62) 



in which we have used Eq. d6*Tl ) and defined s = 9f2 m j42£ 2 - 
Eq. ( f62b is only valid when r ^ 3/2, while the case of r — 
3/2 corresponds to m(a) and /3(a) both being non-power-law, 
which will be studied elsewhere. 

As in the symmetron case, we need to have the expression 
of V v (tp) to study the nonlinear evolution of <p. For this we 
will use the relations 

d(nV) 



da 



= -3 



f3 2 (a) P 2 m (a) 
am 2 (a) 



2s 



2r-3 



(a 



2r-3 



,(63) 



where we have used the expressions of m(a) and /3(a) given 
in Eqs. ( 1611621 ), and 

d(y/Kp) 



= 3 



da 

/3(a) pm(a) 
am 2 (a) M|j 



2 e 2„2r-4 



cxp 



2r 



2r-3 



1 



(64) 



Using the above two equations, it is straightforward to find 

d[nV(a)]/da 



d(y/K(p)/da 
-3fl m /3 HQ exp 



2r - 3 



(a^~ J - 1) 



-6U m n — 



1 , 2r-3 , A 2 (y-y>.) 
s Mpi/3 



a" 3 (65) 



, (66) 



where Eq. ( I65t can be used directly when one needs the back- 
ground value of V v (<p) and Eq. (l66l l can be used in full nonlin- 
ear calculations such as the iV-body simulations. As in general 
A<2.{ip — tp*)/M-p\ < /3q, the logarithmic here is negative, and 
to make sure the last line of Eq. d66l ) is well defined for any r 
we should require r < 3/2. Otherwise the terms in the brack- 
ets can be negative when tp — > tp*, making the power func- 
tion ill-defined. Because tp appears in both (3(p) and V v (p) 
through p — tp*, without loss of generality, in what follows we 
take <y3» = by a redefinition of tp. 



2. Effects of varying model parameters 

As in the symmetron model, let us first analyse how vary- 
ing the four parameters A2,(3o,r and £ affects the structure 
formation. 

1. Increasing enhances s = 90 m A2^ 2 and so makes 
(3(a) smaller at a < 1. As (3(a) controls the strength of 
the fifth force, this weakens its effect. 

2. Increase in (3o makes /3(a) larger at all times, which 
strengthens the fifth force. 

3. The effects of r are two-fold. On the one hand, increas- 
ing r makes m(a) larger and therefore the fifth force 
shorter ranged for a < 1; on the other hand, it makes 
(3(a) larger for a < 1, provided that 2r — 3 is not very 
close to 0, and this strengthens the fifth force. As a re- 
sult, we expect that this will decrease the matter clus- 
tering on large scales but increase it on small scales. 

4. An increase in £ is equivalent to a decrease in mo and 
an increase in s, which means that both m(a) and (3 be- 
come smaller for a < 1. This increases the matter clus- 
tering on large scales and decreases it on small scales. 
Because of the exponential function in (3(a), the effect 
of changing £ is most significant at early times. 

5. There are degeneracies between the different effects. 
For example, increasing r and decreasing £ are expected 
to leave similar imprints on the large-scale structure, as 
we see below. 

Note that the dependence on £ is quite different from that in 
the chameleon models with constant coupling (3 l28ll35[[37ll . 
and the symmetron model lfl7ll . In those cases, increasing £ 
decreases m(a) and therefore increases the range of the fifth 
force, resulting in more matter clustering. 

The above analyses only apply to linear perturbations, the 
dependence of the fifth force on the dilaton parameters is more 
complex in the nonlinear regime, and this is best seen from the 
two functions (3(p) and V v (tp), which govern the nonlinear 
equations (see above): 

1. Increasing Aj implies that the parabolic function A(tp) 
becomes steeper near its minimum at tp = tp*, and this 
makes it harder for the scalar field to roll away from tp*, 
where (3(tp) = 0. This weakens the fifth force. 

2. Increasing (3q makes A2(tp — tp*) / 'Mpi/3o closer to zero 
and therefore |V^,(y)| larger. This means that V(tp) be- 
comes steeper, making it easier for the scalar field to roll 
away from tp* where (3(ip) = and therefore strength- 
ening the fifth force. 

3. If 2r — 3 is not too close to zero, increasing r towards 



3/2 makes 11^,(^)1 larger according to Eq. (1661 ) and so 
makes it easier for the scalar field to roll away from tp* 
where (3(ip) = 0. This strengthens the fifth force. 
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4. Similarly, increasing £ (therefore s) makes V(<p) shal- 
lower and the fifth force weaker. Meanwhile, the scalar 
field becomes less massive and therefore less likely to 
follow the local minimum of V e s which is determined 
by the matter density field and more likely to take larger 
values - this could give rise to a larger value of (3 and 
therefore a stronger fifth force. 



Eq. ( fl~8T > becomes 



(M - N)c 2 S, 2 

„2 



(M - N)5 2 £ 2 



(68) 



IV. THE TV-BODY SIMULATIONS 
A. Equations in code units 

In this section we derive the equations used in the TV-body 
simulations, namely, the Poisson equation for the gravitational 
potential and the EOM governing the dynamics of the scalar 
field. For the sake of completeness we first describe the code 
units used in these equations. The code units used in our code 
are based on (but not exactly the same as) the supercomov- 
ing coordinates of 04511 . They can be summarised as follows 
(tilded quantities are expressed in code units): 



aB' 



a 2 (f> 

(BH y 



pa 

= — , v = 

rr dt . 

-, dt = H ^r, c = 



BH ' 
c 



2. The dilaton case 

Similarly, the dilaton field ip is generally very small (ip <C 
Mpi) and should be positive (otherwise the logarithmic in 
Eq. ( 166b is ill-defined). This means that the numerical value of 
ip can easily go negative in the relaxation procedure, leading to 
the failure of convergence. To avoid this problem, we follow 
||25[ l35ll and use a newly-defined variable u = \og(ip/Mpi) 
instead of ip itself. During the cosmic evolution \u\ remains 
O(l) ~ 0(10), compared to the several orders-of-magnitude 
span of tp, making it easier to handle the numerical errors. 

After some simplification, the dilaton equation of motion 
Eq. (fTST ) becomes 



VV w ^n m A 2 pe u a- 



(69) 



-— n m A 2 e u 
c 2 



2j ,_ 3 2r-3. e" 
1 3 + log — 

s <p 



BH t 







In the above x is the comoving coordinate, p c is the critical 
density today, f2 m the fractional energy density for matter to- 
day, v the particle velocity, <fr the gravitational potential and 
c the speed of light. In addition, B is the size of the simula- 
tion box in unit of /i _1 Mpc and Hq the Hubble expansion rate 
today in units of lOO/i km/s/Mpc. Note that with these con- 
ventions the average matter density is p = 1 at all times. All 
the newly defined quantities are dimensionless. 

Using the code units defined above, the Poisson equation 
Eq. ( fTTI ) becomes 



V 2 0« -n m a(p- 1). 



(67) 



Note that the Poisson equations for both the symmetron and 
the dilaton cases are unchanged compared to the case of stan- 
dard GR, because we have neglected the contribution from the 
scalar field to the source term. In what follows, we introduce 
the symmetron and dilaton versions of the scalar field equa- 
tion, i.e., Eq ( 1181 . 



B. The discretised equations 

Evidently, to put the above equations into the TV -body code 
one must discretise them. For the Poisson equation we have 



Tp \_4>i+l,j,k + <j>i-l,j,k + 0i,j+l,fc + 4>i,j~l,k + <fii,j,k+l 

3 

+^t,j,fc-l - 6< &,j,fc] = 2^™ a (Pi,j,k ~ 1) >(70) 

where 4>i.j.k is the value of <f> in the grid cell with index 

{i,3,k). 



1. Symmetron equation of motion 

The discrete version of the nonlinear symmetron EOM can 
be obtained similarly: 



L h {<p itj ,k) = 0, 
where the operator L h (<pij t k) is defined as 



(71) 



1. The symmetron case 

Throughout the cosmic history, the symmetron field has a 
small magnitude, i.e., \(p\/Mp\ -C 1. To guarantee the numeri- 
cal accuracy, instead of solving ip itself, we solve for a newly- 
defined variable tp = <p/<p*. This variable is constrained by 
< \<p\ < 1 everywhere. The symmetron equation of motion 



L (<Pi,j,k) = T2 [<Pi+l,j,k + <Pi-X,j,k + <Pi,j+l,k + <Pi,j-X,k 
+<fii,j,k+X + ¥>i,j,k-X - 6<Pi,j,k] 



a 



(M - N)d 2 £ 
^2 



,~N-X 



Pi^k-t 



(M - N)5 2 £ 



_~M-1 
7.1t2^i,j,k ' 



(72) 
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Eq. ( TTTb is solved using the nonlinear Gauss-Seidel relaxation, 
which can be summarised as 



~/i,new ~h,old 



where 



(N-l)a 2 N _ 2 



6 

h 2 {M — N)c 2 ^ 2 

(Af-l)q 2 M _ 2 
(M - 7V)5 2 £ 2 ^ ' J '' fc 



In practice, Eqs. (1721741 ) must be modified at the boundaries 
of refinements for the multigrid implementation, as is the case 
of the Poisson equation. Ref. l43ll gives a detailed review of 
all the technical details involved in the TV-body code imple- 
(73) mentation: interested readers are referred to that paper. 



2. Dilaton equation of motion 

The discrete version of the nonlinear dilaton equation can 
be obtained similarly: 



L h {u id , k ) = 0, 



(74) where the operator L h (ui,j tk ) defined as 
I 



(75) 



b iij+ i tk u itj+ i tk - u itj ,h [b i j+ i k + l>, , . , ) + b id _i ik Uij-i,k 



h 2 
1 

'h 2 



c z 



,2r-3 



2r — 3 u 
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+ b ijk _iu itj! k-i 



- —VL m A 2 pi,jMa V 



(76) 



Here b = de u /du 



bi+±,j,k = 9 {h+i,j,k + bi,j,k) 
b 



and /i is the length of the cell in the numerical simulation 
mesh. 

I 



I 

Eq. (l75T l is solved using the nonlinear Gauss-Seidel relax- 
ation as well, which can be summarised as 



/i,new /i.old 
U, '■ u = U. ' 



i,j,k "i,j,k dLl(ulf£) ' 
()u 1 



,h, old 



where 



(77) 



dL h (uij.k) 



du 



-^Tgbi^fr [lij+i^fc + Ui-ij-,fe + Ui,j+i,/c + Ui-\.j ik + Ui dt k+i + u i,j,k-i — 6ui,j,fc] 

e 2 r , 



-3£l m j4 2 e 



2r — 3 it 



2„«; 



-3£l m ^4 2/ 5a e 



2r-3 , 2r 3 Ujj^fc 



(78) 



Again, Eqs. d76l l and (T78b must be modified at the bound- 
aries of refinements for the multigrid implementation, as is the 
case of the Poisson equation. 



V. CODE TESTS 

In this section we present the results of code tests we have 
performed to show that our symmetron and dilaton equation 
solvers work well. To lighten the notation, throughout this sec- 
tion we use the units AIpi = 1. 

There are five parameters for the generalised symmetron 
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TABLE I. The parameter values for the six models used in the sym- 
metron code test. 



model a* A) (N, M) 

a 05 05 (274) MOOT 

b 0.2 0.5 (2,4) 0.001 

c 0.5 1.0 (2,4) 0.001 

d 0.5 0.5 (2,6) 0.001 

e 0.5 0.5 (2,4) 0.0005 

f 0.5 0.5 (M) 0.002 



TABLE II. The parameter values for the five models used in the dila- 
ton code test. 

model A2 /3o r f; 

a 5 x 10 b 05 I 0.001 

b 1 x 10 6 0.5 1 0.001 

c 5 x 10 5 1.0 1 0.001 

d 5 x 10 5 0.5 0.001 

e 5 x 10 s 0.5 1 0.002 



model, namely a*, f3 , N, M and £, and we set N = 2 and test 
the code for 6 models summarised in table U There are 4 pa- 
rameters for the generalised dilaton model, namely A2, (3q, r 
and £ (note that s can be calculated when A2 and £ are given, 
and is therefore not an independent model parameter), and we 
test the code for 5 models summarised in table IIT1 



s- 

1E-5 





I 


1 


I 1 


1 1 1 
















- 


model a. 


a=1.0, 


before rel. 


model a, a~1.0, 


■ 

after re/. 




model b, 


a=1.0 


before rel. 


model b, a=1.0, 


after re/. 




model a, 


a=0.6 


before rel. 


model a, a=0.6, 


after re/. 




□ □ □□ 






==" nn° D on* ° "= 




XA □ 






□ □□ 

M * 
........... 








I 




1 


A A 
I , I 





0.0 0.2 0.4 0.6 0.8 1.0 



X/B 



FIG. 3. (Colour online) Test of the solver for the symmetron equation 
in a constant matter density field. Only results in the cells along the 
a;-axis are shown, and the x-coordinate is rescaled by the size of 
the simulation box so that x G [0, 1]. Results for three models as 
explained in the legend have been shown (the empty symbols), the 
final answer corresponding to which are filled symbols of the same 
type and colour. The horizontal lines with the same colours are the 
exact analytical solution. 



A. Homogeneous matter density field 

In a universe with a homogeneous density, the symmetron 
field ip should exactly take its background value ip, namely 



<f {a) = </?* 



(79) 



everywhere. Thus, as the simplest test of the symmetron equa- 
tion solver, one can show that in such a homogeneous field, 
given some random initial guess of <p on the cells of the sim- 
ulation mesh, after a reasonable number of Gauss-Seidel re- 
laxation sweeps, the solutions all converge to the above back- 
ground value. Such simple test have been used previously in 
ll39l I4TL l43ll to show that the solver for extra degrees of free- 
dom works correctly. 

We have performed this test for all the six symmetron mod- 
els summarised in TableU The result is shown in Fig. [5] where 
we plot the values of ip/Mp\ in the cells in the x-direction, be- 
fore and after the Gauss-Seidel relaxation; for clarity we have 
only shown the results for models a and b at a = 1.0 and 
model a at a = 0.6. We can see that the final solution agrees 
with the analytical result (the horizontal lines) very well (see 
figure caption for more details). 

We have also tested the code for a model with a* = 0.5 at 
a = 0.4. In this case the symmetry of V c s((p) has not been 
broken yet, and we expect that ip vanishes everywhere. This is 
confirmed by the tests (which are not shown here). 



For the dilaton model, the field ip also takes exactly its back- 
ground value ip, given by 



— / \ A) - 

A.2 



exp 



,2r-3 



2r - 3 



(80) 



everywhere in a homogeneous universe. 

We have performed this test for three of the five models 
summarised in TablellTl The results are shown in Fig. [4] where 
we plot the values of log(ip/(p) in the cells in the x-direction, 
both before and after the relaxation. For clarity we have only 
shown the results at a — 1.0. It can be seen that the final solu- 
tion agrees with the analytical result (the horizontal lines) very 
well (see figure caption for more details). We have also tested 
our code at a 7^ 1.0 and found the same good agreement. 



B. Point mass 

As a second test of our symmetron equation solver, let us 
consider the solution of ip around a point mass at the origin, 
for which case we have an analytical solution which is accu- 
rate except for the regions very close to the mass. Such a test 
has been used previously in ll25l . l39l.l43l] . 

Following [25], we construct the point-mass density field as 
(hereafter S idtk = p i>jik - 1) 



10~ 4 (N 3 -1) , i = j = k = 0; 
— 10~ 4 , otherwise. 



(81) 
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FIG. 4. (Colour online) Similar to Fig. [3] but for the dilaton model. 
For clarity only the results of models a, b, d (as indicated in the leg- 
end) are shown: the initial guesses are represented by the empty sym- 
bols and the numerical solutions are denoted by filled symbols of the 
same type and colour. Note that, instead of \og(ip), we have shown 
\og(ip/if). The horizontal lines with the same colours are the exact 
analytical solution, which is zero identically. 



in which i, j, k are respectively the cell indices in the x, y, z 
direction. In the test we use a cubic box with size 250/i~ 1 Mpc 
and 256 grid cells in each direction. We have done this test for 
all six models of table U at a = 1. 

On the other hand, the analytical solution can be obtained 
approximately by solving the equation 



m 2 Sip 



(82) 



in which the effective mass of the scalar field Sip — <p — (p is 



! tt2 

H . 



The analytical solution is 

Sip oc — exp(— mr), 
r 



(83) 



with r the distance from the point mass. 

Fig. |5] shows the comparison between the numerical solu- 
tions to S(p along the a>axis (symbols) and analytical solutions 
(solid curves) for the symmetron models, and we can see that 
the two agree very well in all cases. The discrepancies at small 
x is because the linearisation procedure in deriving Eq. ([821 is 
not accurate and the discrepancy at big x is because the size of 
Sip has reached the level of the discretisation error [|25l . Fig. [6] 
shows the comparison for the dilaton models, and once again 
we find excellent agreements. 



C. Sine density field 

As our third test, let us consider the sine density field intro- 
duced in i25ll . which (after some modification to account for 



FIG. 5. (Colour online) The solution to Sip = ip — ip around a point 
mass constructed according to Eq. (18 H . for the six test symmetron 
models in Table U (see the legend). The solid curves with the same 
colours are the corresponding analytical approximations which are 
accurate far from the point mass. Only solutions along the x-axis are 
shown. 



the code units) in the symmetron case is given by 



(M - N) sin(27ra) 



[2 - sin(27rx)] JV - 
- [2 -sin(2 7 ra;)] M - Ar , 



(84) 



where x is rescaled so that x <G [0, 1]. We consider only the 
x-dependence, which is equivalent to a one-dimensional con- 
figuration. The solution to this density field can be analytically 
worked out to b^| 

ip{x) = tp*[2 - sin(27ra)]. (85) 

Fig.|7]shows the symmetron test results for the sine density 
field given above, at a = 1 and for the six models listed in 
Table U It can be seen that the numerical solutions (symbols) 
agree with the analytical solutions (solid curves) very well. 

Similarly, for the dilaton field let us consider the following 
density field 

c 2 a (2tt) 2 sin(27rx) 



p(x) 



n m A2 3 2 — sin(27ra;) 



(86) 



,2r-3 



2r - 3 



loe 



sin(27ra;) 



in which x is rescaled such that x 6 [0,1]. The solution to this 
density field can be analytically worked out to be, 



ip(x) 



&2 



i(2irx)} . 



(87) 



6 More exactly speaking, we specify the solution we want the code to repro- 
duce and then use the EOM to calculate the corresponding density field that 
gives rise to this solution. 
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FIG. 6. (Colour online) The solution to Sip = cp — <p around a 
point mass constructed according to Eq. (18 H . for the five test dila- 
ton models in Table ITT1 (see the legend). The solid curves with the 
same colours are the corresponding analytical approximations which 
are accurate far from the point mass. Only solutions along the rc-axis 
are shown. 

Fig.|8]shows the dilaton test results for the sine density field 
given above, at a = 1.0 for models a, b, c and at a = 0.2 
for model a listed in Table [II] As in the symmetron case, the 
agreement is very good. 




0.0 0.2 0.4 0.6 0.8 1.0 



X/B 



FIG. 7. (Colour online) Solutions of ip in a one-dimensional (re- 
direction) sine density field constructed using Eq. (HD, for the six 
test symmetron models (as indicated besides the curves). The solid 
curves with same colour are the corresponding analytical results and 
the symbols are the numerical solutions. A simulation box with side 
length of 250ft -1 Mpc and 256 grid cells on each side is used in the 
computation, x is rescaled so that x/B G [0, 1]. 



D. Gaussian density field 

The last test on the regular (i.e., unrefined) grid uses a Gaus- 
sian type density configuration. Again, here we only consider 
one dimension, and for the symmetron case the density field 
is specified as 



(^)'p(*) = l + (^ 
V a J \ a 



1 — a exp 



a(M - N)(x - 0.5) 2 /W 2 



^1 — a exp 
(x-0.5) 2 



(x-0.5) 2 
W 2 



M-N 



W 2 



N-l 



(88) 



where again x has been scaled to code units so that x E [0, 1], 
W, a are numerical constants which respectively specify the 
width and height of the density field, which obviously peaks 
at x = 0.5. Such a density field has been used in the code test 

of m. 

Note that such a density field is not exactly periodic at the 
edges of the simulation box, but given that W is small enough, 
p —> at the box edges and periodic boundary conditions are 
approximately satisfied. 

The solution to cp can then be obtained analytically and is 



ip(x) = ip* 



1 — a exp — 



(x-0.5) J 
W 2 



(89) 




10- 7 , 
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model a, a=0.2 
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FIG. 8. (Colour online) Solutions of ip in a one-dimensional (re- 
direction) sine density field constructed using Eq. d86) . for three test 
dilaton models (a, b, c) at a — 1.0 and model a at a = 0.2 (as indi- 
cated besides the curves). The solid curves are the corresponding an- 
alytical results and the symbols are the numerical solutions. A simu- 
lation box with side length of 250/i~ 1 Mpc and 256 grid cells on each 
side is used in the computation, x is rescaled so that x/B 6 [0, 1]. 
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FIG. 9. (Colour online) Solutions of <p in a one-dimensional (re- 
direction) Gaussian-type density field constructed using Eq. J88b . for 
the six test symmetron models in Table U (see legends). The solid 
curves are the analytical results from Eq. (1891 and the symbols with 
same colours are the corresponding numerical solutions. A simula- 
tion box with side length of 250/i _1 Mpc and 256 grid cells on each 
side is used in the computation and the symmetron equation is only 
solved on the regular domain grid, x is rescaled so that x/B £ [0,1]. 



FIG. 10. (Colour online) Solutions of <p in a one-dimensional (re- 
direction) Gaussian-type density field constructed using Eq. ( 1901 ), for 
three test dilaton models (a, b, c) at a = 1.0 and test model a at a = 
0.3 (see legends). The solid curves are the analytical predictions from 
Eq. ( 1891 ) and the symbols with same colours are the corresponding 
numerical solutions. Other specifications are the same as in Fig. [9] 



which clearly shows that when a 1 \ip\ could be made very 
small at x = 0.5 while at x — > or x — > 1 it goes to cp = (p*. 

We have implemented Eq. d88l ) into our numerical code and 
the numerical solutions for ip are shown in Fig. [9] We can see 
that they agree with the analytical solution Eq. d89b very well. 

For the dilaton case we use the following density field 



p{x) 



~c 2 a 2a ex P 



(x-0.5) 2 
W 1 



sn m A 2 w 2 



1 — a exp 

2r-3 



•i - * — w 1 — 



(z-0.5) 2 
W 2 



(90) 



■log 



1 — ae 



3 

2r-3 



where x, W and a are specified similarly as above. 

The test results for the dilaton models are shown in Fig.fTUl 
where again we find good agreement with the analytical solu- 
tion Eq. 



E. Equation solver on refinements 

The above tests show that our solver of the scalar field EOM 
works accurately on regular grids. But in cosmological simu- 
lations these equations are also solved on irregularly-shaped 
refinements where they can take different forms due to the re- 
finement boundaries H43fl - It is therefore necessary to test the 
scalar field equation solver on refinements as well, which we 
will do in this subsection. 




a=0.999, level 8 
a=0.999, level 9 
a=0.9999X0.1, level 
a=0.9999X0.1, level 9 
a=0.99999X0.01, level 8 
0=0.99999X0.01, level 9 
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FIG. 11. (Colour online) Same as Fig. [5] but for the model a only and 
a = 0.999, 0.9999, 0.99999 (from top to bottom: red, green, blue). 
The symmetron equation is solved on two levels: level 8 (the regular 
domain grid) and level 9 (the first refinement), and their numerical 
solutions are represented by empty and filled symbols of the same 
shape and colour respectively. The solid curves of the same colours 
are the corresponding analytical solutions from Eq. ( 1891 . A simula- 
tion box with side length of 250/i~ 1 Mpc and 256 grid cells on each 
side is used in the computation and the symmetron equation is only 
solved on the regular domain grid, x is rescaled so that x/B £ [0, 1]. 
For clarity we have multiplied the results for a — 0.9999 and 
0.99999 by 0.1 and 0.01 respectively. 
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The Gaussian-type density configuration provides a good 
way to check the multilevel scalar-equation solver, because 
the density peak can be made arbitrarily high by adjusting the 
parameter a and the value of the matter density is the crite- 
rion we use to refine grid cells in cosmological simulations. 
In the vicinity of this peak, the density field p changes rapidly 
and higher spatial resolution is necessary to compute ip (and 
differentiate it to get the fifth force) accurately. 

Consider the case where the regular domain grid is refined 
only once, in regions where the density value exceeds a given 
threshold (we call this a 'two-level problem', and in the nu- 
merical examples below the coarse and fine levels are respec- 
tively levels 8 and 9). The density values p in both the coarse 
and the refined cells are calculated using Eq. (l88l for the sym- 
metron case and Eq. d9Qb for the dilaton case, while the values 
of ip at the fine-level boundaries are computed from interpola- 
tion of those in the nearby coarse-level cells ll43ll . 

Fig. [TTI shows the numerical values of ip on both levels in 
the region covered by the refinement, for the symmetron case. 
We show the results for model a only and for four different 
values of a (0.999, 0.9999 and 0.99999 from top to bottom), 
and for each a the results from the coarse and fine levels are 
denoted respectively by empty and filled symbols. For com- 
parison we have also plotted the analytical results Eq. (l89l as 
solid curves. As we can see, both fine-level and coarse-level 
results are virtually indistinguishable from the exact solution. 

This does not mean that the refinement is unnecessary how- 
ever, because, as shown in Fig.QT] the fine level has more data 
points and could probe regions closer to the extreme value of 
(p, which corresponds to the high density region where high 
resolution is needed. 

For the dilaton, Fig.Q~2]shows the numerical values of (p on 
both levels in the region covered by the refinement. Again, we 
show the results for model a only and for four different values 
of a (0.999, 0.9999 and 0.99999 from top to bottom), and for 
each a the results from the coarse and fine levels are denoted 
respectively by empty and filled symbols. For comparison we 
have also plotted the analytical results Eq. (l89b as solid curves. 
Excellent agreement is found again. 



F. Other tests 

In the above we have focused on various tests of the scalar 
field solver of the ECOSMOG code, as this is the only new addi- 
tion to the default RAMSES A-body code. These tests checked 
the validity of the new subroutines against different density 
distributions, and the good agreements with analytical solu- 
tions shows the validity of the code and its accuracy. 

As the standard gravity solver and particle-updating sub- 
routines of RAMSES are not touched, tests carried out for them 
(which show that the RAMSES code works very well) need not 
be repeated here. The AMR simulation algorithm is often im- 
plemented in different ways in different codes; for a detailed 
explanation of its implementation in RAMSES and therefore in 
ECOSMOG we refer to Hi and |0] respectively. We do not 
present the full details here as they are too long and this paper 
is mainly concerned with the modified gravity physics. 
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FIG. 12. (Colour online) Same as Fig.fTO] but for the model a only 
and q = 0.999, 0.9999, 0.99999 (from top to bottom: red, green, 
blue). The dilaton equation is solved on two levels: level 8 (the regu- 
lar domain grid) and level 9 (the first refinement), and their numerical 
solutions are represented by empty and filled symbols of the same 
shape and colour respectively. The solid curves of the same colours 
are the corresponding analytical solutions from Eq. ( 1891 . A simula- 
tion box with side length of 250/i _1 Mpc and 256 grid cells on each 
side is used in the computation and the dilaton equation is only solved 
on the regular domain grid, x is rescaled so that x/B £ [0, 1]. For 
clarity we have multiplied the results for a — 0.9999 and 0.99999 
by 0.5 and 0.25 respectively. 

When a new code is written, one needs to test its cosmolog- 
ical simulations. This is straightforward for a standard code 
of ACDM simulations, because there are fitting formulae and 
results from other codes to compare to. Unfortunately, up to 
now there are no accurate fitting formulae for modified grav- 
ity theories such as symmetron, dilaton and f(R) gravity. But 
several serial A-body codes simulating f(R) gravity (e.g., 
||25[ HU) and symmetron models (e.g., Mill ) do exist in the 
literature: in both cases good agreement with ECOSMOG has 
been founcfl See, for example, [43] for a comparison for f(R) 
gravity, and we have also checked explicitly that our sym- 
metron simulation result agrees with that of 114 ill . 

Finally, for cases where approximate analytical results can 
be obtained from other methods, we find good agreement be- 
tween ECOSMOG and the approximation solutions. An exam- 
ple is the f(R) gravity model of 10 with |d//dR| = 10~ 4 , 
the nonlinearity of which is very weak and so the matter power 
spectrum can be approximated by linear perturbation theory 
down to relatively small scales. This is actually confirmed in 
ll34ll . which can serve as another test of the ECOSMOG code. 

In short, the ECOSMOG scalar field solver has been tested in 



7 Another independent code which is still being developed also agrees with 
ECOSMOG very well. 
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various ways, and several cosmological simulations of modi- 
fied gravity models using ECO SMOG agree with similar sim- 
ulations done using other codes, such as the codes developed 
independently in HI SH . 



VI. COSMOLOGICAL SIMULATIONS 

In this section we describe and analyse the results of cos- 
mological simulations of the dilaton and symmetron modified 
gravity models. We also perform ACDM simulations for com- 
parison. For each model we run 5 realisations with the same 
physical parameters and simulation specification, but differ- 
ent realisations of initial conditions. The initial conditions are 
generated using MP GRAF I C JH] at redshift Zi = 49.0 with 
different seeds of random numbers. Since at z, = 49.0 the ef- 
fect of the fifth force is negligible, the initial conditions should 
be the same for all models studied here. For the ease of com- 
parison, we use the same random seed to generate initial con- 
ditions for the same realisation of all models, including sym- 
metron, dilaton and ACDM. 

The background expansion history in the studied dilaton 
and symmetron models is in practice indistinguishable from 
that of the fiducial ACDM model I17I1 . In all simulations we 
adopt WMAP7 14711 cosmological parameters, with h = 0.71, 
n m = 0.267, n A = 0.733, n s = 0.963 and cr 8 = 0.801. 

The size of the simulation box is chosen to be 128/i _1 Mpc, 
and the domain gricflhas 2 8 = 256 cells on each side. The grid 
cells are refine when the effective number of particles in them 
exceeds 9.0, and the finest refinement level equivalently has 
2 14 cells on each side. The number of particles is N p = 256 3 
in all simulations. 



A. The symmetron models 

The symmetron models are specified by the four model pa- 
rameters a*, M, N and £. We have chosen to fix /3* = 1.0 
for all our runs in order to see the effect of varying the other 
parameters individually. The effect of varying /3+ is to modu- 
late the strength of the fifth force and was investigated for the 
symmetron in J4lll . In Table (ITXTb we list the parameters for the 
nine models we have simulated. 

In the rest of this subsection, we will focus on the effects 
of changing each model parameter on the major cosmological 
observables such as the matter power spectrum and halo mass 
function. More specifically, we will analyse the results of our 
numerical simulations according to the following: 

1 . How the symmetry breaking scale factor a* affects the 
results: Model Al versus Bl, A2 versus B2 and A4 ver- 
sus B4. 



TABLE III. The parameter values for the nine models used in the 
symmetron cosmological simulations. For each model we have 5 re- 
alisations of initial conditions, and therefore a total of 45 runs. 



model name 


a* 




(N, M) 


2998£ 


realisations 


ACDM 


— 






— 


5 


Al 


0.50 


1.0 


(2,4) 


1.0 


5 


A2 


0.50 


1.0 


(2,6) 


1.0 


5 


A3 


0.50 


1.0 


(2,6) 


2.0 


5 


A4 


0.50 


1.0 


(4,6) 


2.0 


5 


Bl 


0.33 


1.0 


(2,4) 


1.0 


5 


B2 


0.33 


1.0 


(2,6) 


1.0 


5 


B3 


0.33 


1.0 


(2,4) 


2.0 


5 


B4 


0.33 


1.0 


(4,6) 


2.0 


5 



2. How the coupling strength parameter N affects the re- 
sults: Model A3 versus A4 and B3 versus B4. 

3. How the potential parameter M influences the results: 
Model Al versus A2 and Bl versus B2. 

4. How the range A* = 2998£ Mpc /h of the fifth force in- 
fluences the results: Model A2 versus A3 and Bl versus 
B3. 



1. Nonlinear matter power spectra 

The most direct way to see the effect of modified gravity 
on the clustering of matter is to look at the matter power spec- 
trum P(k). We have measured the nonlinear P(k) in the sym- 
metron models and calculated their relative differences from 
the ACDM prediction. The results are shown in Figs. [T3lfT4l 
The power spectra are measured using the publicly available 
code POWMES JH]. 

1. The symmetry breaking scale factor a* controls when 
the fifth force starts to kick in. From Fig. Qj] we could 
see that decreasing a* (i.e., moving from A models to 
B models) leads to a stronger matter power spectrum as 
the fifth force would have more time to participate in 
structure formation. Notice that when a < a* the mat- 
ter power spectra in symmetron models are essentially 
unchanged as can be seen in Fig. [T-ffl . This is because 
on linear scales there is strictly no fifth-force effect be- 
fore a = a*, since the magnitude of the fifth force is 
determined by the background matter density, which 
is always higher than p* before a = a*. However, on 
nonlinear scales, the fifth force can kick in even before 
a = a* in regions where matter density drops below p*, 
thus the structure formation is affected even at a*. 



As RAMSES and ECOSMOG are adaptive mesh refinements codes, the do- 
main grid is defined as the finest uniform (regular) grid which covers the 
whole simulation box. 



9 In Fig. ll4l svmmetrv breaking has just happened for A models and the fifth- 
force effect has not accumulated at a = 0.5. 
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FIG. 13. (Colour online) The relative difference between the matter power spectra of the symmetron models and the ACDM paradigm. The 
symbols are from the iV-body simulations, and the curves are linear perturbation theory predictions. Details are illustrated by the legends, and 
a = 1.0. 



2. The parameter N of the matter coupling j3 oc t^A" 1 de- 
termines how the matter coupling evolves. As the field 
moves towards ip = in high density regions, a larger 
N means that the fifth force becomes more suppressed 
as shown in fl7il . This effect can be seen in Fig.ll3l(up- 
per right panel). Note that varying N also changes the 
evolution of ip through the changes of (3(<p) and V(<p); 
however the numerical result here shows that this effect 
is subdominant. 

3. The parameter M of the self-interaction term <p M E 
V(ip) determines how nonlinearly the model behaves. A 
higher-order (larger M) interaction term means that the 
nonlinearities, and therefore the screening mechanism, 
are less at play (see § 1111 A 2b . which again leads to more 
matter clustering as confirmed by the lower-left panel of 
Fig. [T3] This effect can also be seen by noting that the 
nonlinear power spectra for the cases of M = 6 are in 
general closer to the corresponding linear power spectra 
than for the cases of M = 4. 

4. The range A* = 2998£Mpc//i of the fifth force deter- 
mines which scales are influenced by the fifth force. In- 



creasing the range moves the modifications of gravity 
to larger cosmological scales as can be seen in Fig. [13] 
In the linear perturbation regime, the power spectra 
for two models with different ranges (A*]^) are re- 
lated by the scaling relation Pi(k) — PiikA+i/ X^)- 
However, this scaling no longer holds in the nonlin- 
ear regime. For example, when A* decreases, the sym- 
metron mass becomes heavier, the screening effect is 
enhanced and consequently the power spectrum is sup- 
pressed (c.f. Fig.[T3]and § 1111 A 2b . 

5. At late times (Fig.Qj) the linear perturbation prediction 
is a bad approximation to the full solution, which is be- 
cause the symmetron EOM is highly nonlinear. Indeed, 
as in the case of f(R) gravity [34], the linear theory be- 
comes inaccurate almost as soon as the power spectrum 
starts to deviate from the ACDM prediction. This shows 
the importance of properly taking into account the non- 
linear effects (by numerical simulations) in the study of 
structure formation in modified gravity models. 

6. The agreement between the linear and nonlinear results 
becomes better at earlier times (c.f. Fig. [14), when the 
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effect of nonlinearity has not accumulated for long. 

In f(R) gravity models, it is known l34ll that the shape 
of AP/P follows a fixed evolution path, and at any given 
time the position of a model on this path is determined by the 
properties of the fifth force and how long it has become non- 
negligible. Similar patterns appear here, for example in the A 
models a* = 0.5 where the fifth force becomes non-negligible 
later than in the B models, for which a* = 0.3. Correspond- 
ingly, in Fig. Qj] AP/P has a peak at k — l/iMpc -1 . On the 
other hand, Fig.[T3lshows that for symmetron models AP/P 
goes up again on very small scales (k > a few), while in f(R) 
models AP/P decreases for these scales 0411 . 

This pattern for the symmetron matter power spectrum can 
be understood as follows. At early times the model is well 
described by linear perturbation theory and the symmetron 
mass (and the coupling strength f3(ip)) is nearly the same ev- 
erywhere; the Yukawa nature necessarily means that the fifth 
force decays with distance, and as a result AP/P increases 
monotonically with k at these times (see Fig.[T4b. Later when 
highly nonlinear and dense structures have formed, the sym- 
metron screening mechanism starts to work so that the fifth 
force inside these structures are efficiently suppressed (/3(<p) 
becomes small) and GR is locally restored since then, which 



makes AP/P frozen on small scales (thus remain monoton- 
ically increasing) while at the same time still grow on larger 
scales (e.g., k > 1/iMpc" 1 ) as the fifth force still propagates 
among different halos. 

To understand this behaviour more properly would require a 
detailed study of the density and velocity fields, together with 
their time evolutions, and these will be left to future work with 
higher-resolution and larger simulations. 

As an illustration of the above effects, the difference be- 
tween the symmetron models we have simulated and ACDM 
on scales of order 1 Mpc can be as large as 30 percent today. 
This can be seen in Fig. 13 for models Bl and B3 where the 
range of the force is respectively 1 and 2 Mpc and the highest 
power in the potential is 6 and 4 respectively. On these ex- 
amples, the characteristic bump of the symmetron models can 
also be seen in a clear way. 



2. Mass functions 

We have measured the mass functions from our simula- 
tions using the publicly available code AHF l49tl . which is ef- 
ficiently parallelised using MP I and OpenMP. The mass of a 
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FIG. 15. (Colour online) The ratio between the mass functions of the symmetron models and the ACDM paradigm at a = 1.0. 



halo is defined as the total mass contained in i?2oo> the radius 
at which the density contrast A drops below 200 times the crit- 
ical density. For each model, including ACDM, we have cal- 
culated the average and standard deviation of the mass func- 
tion over the five realisations. 

Because we are interested in how the fifth force can change 
the matter clustering, we show the ratio of the symmetron and 
ACDM mass functions, 1Z = n symmc tron/^ACDM- The stan- 
dard deviation an of 1Z for each mass bin is computed using 
the normal rule of propagation of errors, according to which 
we have 
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The subscripts mg and a denote the modified gravity model 
(the symmetron here and the dilation in the next section) and 
ACDM respectively, and p is the correlation coefficient be- 
tween the mass functions of the two, i.e., 
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where the sum is over five realisations and the quantity with 
an overbar denotes the average over five realisations. 

In Fig.[T5]we show the ratios between the symmetron and 
ACDM mass functions from our simulations at a = 1.0. The 
results at a = 0.5 are shown in Fig. [16] 

The fifth force leads to an overall enhancement of the for- 
mation of dark matter structures. The effect is strongest for 
intermedium-sized (M ~ 10 13 /i _1 Af Q ) halos and we find a 
maximum enhancement in the mass function of around 50% 
compared to ACDM for the models we have simulated. For 
the largest halo masses (M > 1O 14 /i _1 M ) the symmetron 
mass function goes towards ACDM as the symmetron screen- 
ing mechanism makes sure the fifth force is effectively sup- 
pressed for such massive objects. 

The effects of varying different model parameters on the 
mass function are not as clear as in the power spectrum, but 
we can see the same trends. More specifically, 

1. For models with smaller a* (i.e., the B models) we see 
from Fig. [T5]that a larger fraction of high mass halos is 
obtained. As with the matter power spectrum, the mass 
function is essentially unmodified for a < a* (see A 
models in Fig. [16] for which a* = 0.5 and the effect of 
the fifth force has not accumulated at a = 0.5). These 
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FIG. 16. (Colour online) The same as Fig.[l5]butfora = 0.5. 



are to be expected since the fifth force is not at play on 
cosmological scales at such early times, and for smaller 
a* the fifth force has acted for a longer period. Hence 
more large halos form and fewer small halo survive the 
mergers. 

2. As mentioned in § IIII A 21 increasing the parameter N 
leads to a suppression of the fifth force, especially for 
large halos and in high density regions where \<p\ <C 
This can be seen from the upper-right panel of Fig. Q3] 
Note that in models B3 and B4 both N and M are dif- 
ferent, and the effect is not purely due to varying N. 

3. As discussed in § 1111 A 21 increasing M makes it easier 
for the scalar field to roll away from ip = where the 
coupling strength vanishes. This leads to a stronger fifth 
force and consequently more large halos, as can be seen 
in Fig. IT61 (lower-left panel). 

4. Increasing £ increases the range A* of the fifth force 
and leads to more high-mass halos. This can be seen in 

Kg.na 

As for AP/P, the effects of varying different model pa- 
rameters on the shape of An/n are similar, which shows that 



the four parameters are highly degenerate. This behaviour is 
different from what we will see in the dilaton simulations be- 
low. 

The significant deviations of our symmetron models from 
the prediction of the ACDM paradigm, as shown in Figs. [T51 
andQ~6] should be detectable by future surveys. 



B. The dilaton models 

In this subsection we analyse cosmological simulations of 
the generalised dilaton models. We vary all four model param- 
eters A2, (3o,r and f, so that each of them takes 4 (3 for Ai) 
different values with the rest remaining the same. This results 
in a total of 12 dilaton models, as summarised in Table ITVl 
The choices of parameter values are such that A2, B2, C2 and 
D2 are the same model, to facilitate a cross comparison. 

As the dilaton simulations were run on a different machine 
from the symmetron ones, we have simulated the same default 
ACDM models on both machines, and checked that they agree 
very well. This enables a direct comparison between dilaton 
and symmetron simulations if needed. 
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FIG. 17. (Colour online) The relative difference between the matter power spectra of the dilaton models and the ACDM paradigm. The symbols 
are from the TV-body simulations, and the curves are linear perturbation theory predictions. Details are illustrated by the legends, and a = 1.0. 



TABLE IV. The parameter values for the 65 cosmological simula- 
tions we have performed for this study. Note that '-' means that the 
parameters are unused for the ACDM case, and it means that the 
parameters are the same as in A2 in the cases of B2, C2 and D2. 

model name A2 fto r £ realisations 

5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 



1. Nonlinear matter power spectra 

This subsection contains results about the nonlinear matter 
power spectra for the simulated dilaton models. Fig.fTTlshows 
the relative differences between the dilaton and ACDM results 
at a = 1.0, from which we can see the following properties: 

1 . Decreasing A2 leads to stronger matter clustering, since 
A-2 controls the steepness of the coupling function A(<p) 
(see Fig. [TJ. As discussed in § 1111 B 21 the larger A2 be- 
comes, the steeper A(<p) is and the harder it is for ip to 
roll away from tp* where f3(ip) = - this means that 
f3 is closer to zero and the fifth force is more strongly 
suppressed. 

2. Increasing (3q leads to stronger matter clustering, as /?o 
determines the strength of the fifth force. 

3. The r-dependence is weak since large changes in f3 only 
take place at early times (see below). We see the feature 
discussed in § IHIB 21 that increasing r decreases the 
matter power on larger scales (k < 0.2Mpc/h) and in- 
creases it on smaller scales; this happens in both linear 
and nonlinear results. 
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FIG. 18. (Colour online) The same as Fig.[l7]butfora = 0.5. 



4. As discussed in § IIII B 21 decreasing £ simultaneously 
increases the strength and decreases the range of the 
fifth force, causing more (less) clustering of matter on 
small (large) scales. This can be seen by comparing the 
results of Dl and D2. On even smaller scales, however, 
the matter power spectrum increases with £ again. 

5. As in the symmetron case, at late times the linear per- 
turbation theory is a rather bad approximation to the full 
nonlinear dilaton model, and fails to accurately predict 
the matter power spectrum even for k ~ 0.04/i/Mpc. 
This once again shows the important role iV-body sim- 
ulations have to play in the studies of modified gravity 
theories. 



f(R) gravity models), for which only the first part con- 
tributes to the screening. 

At a = 0.5 (cf. Fig. [18), all the above properties remain, 
with the following noticeable features: 

1. The agreement between linear perturbation theory and 
the full simulations gets better as nonlinearities have not 
reached their full effect. This is the same as the sym- 
metron (see above) and f(R) 13411 cases. 

2. The difference between the different C models becomes 
larger than at a = 1.0 because, as mentioned above, the 
effect of changing r is mainly to modify [3(a) at early 
times. 



6. Overall, we see that the nonlinearity suppresses the mat- 
ter power compared with the linear theory predictions, 
which shows that the dilaton mechanism works well for 
large scale structures. The suppression of the fifth force 
comes from two parts: the smallness of ip and therefore 
V<y9 in high density regions, and the smallness of f3(<p) - 
this indicates that with the same configuration of ip the 
fifth force in the dilaton models here is more strongly 
suppressed than in the case of a constant f3(ip) (e.g., in 



The linear-nonlinear agreement is even better at a = 0.3 
(see Fig . |T9j . This indicates that the nonlinearity of the model 
only becomes important at late times, which is possibly be- 
cause the formation of high density structures only then drives 
(p to deviate from its background value. 

Most of our simulation results show less deviation between 
the simulated dilaton models and ACDM than the case of the 
symmetron models. One of the reasons for this lies in the sim- 
ulation details. In the symmetron models we have fixed the 
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FIG. 19. (Colour online) The same as Fig.[l7]butfora = 0.3. 



coupling strength /?* = 1, while for the dilaton cases, except 
for models B3 and B4, the coupling strength is taken to be at 
most flo < 0.5. As the fifth force scales as /3 2 , this makes a 
significant difference (c.f. Fig. [17] upper right panel). As an 
example, model B4 differs from ACDM by nearly as much as 
the symmetron models do (and even more). 

The shapes of the dilaton matter power spectra are worth 
discussing, as they show significant difference from the cases 
of symmetron and f(R) gravity models. From Figs. [P71 [T8l 
and[19]we can see that: 

1 . In both linear and nonlinear cases, AP/ P tends to flat- 
ten on small scales. In the linear case, this is very differ- 
ent from the behaviour of chameleon models with con- 
stant coupling strength (3. In that case, the fifth force al- 
ways has the same strength but at early times its range is 
limited by the very heavy scalar field mass: this means 
that on very small scales the fifth force has started en- 
hancing clustering of matter ever since very early times, 
which is why AP/ P keeps increasing with k jlln . For 
dilaton models, on the other hand, the scalar field mass 
evolves more slowly and the coupling strength is sup- 
pressed at early times: this means that by the time the 
fifth force becomes non-negligible, its range has be- 



come large enough and below this range the growth 
of matter density perturbations is enhanced in a nearly 
scale-independent way (at least in the linear regime). 
Such a feature can indeed also be seen in the linear pre- 
dictions of AP/ P for symmetron models (cf. Fig. [13). 

2. The flattening effect of AP/P on small scales is pre- 
served when varying model parameters A2 and /3 U , but 
is weakened by varying r and £. This is because, as dis- 
cussed in § IIII B 21 varying A 2 and /3 does not change 
the scalar field mass to, while varying the other two pa- 
rameters does. Taking the parameter r as an example, 
increasing r makes m more sensitively dependent on 
local matter density (i.e., more like a chameleon model 
which has no flattening in AP/P). On the other hand, 
decreasing r makes j3 more sensitively dependent on 
local matter density and so suppresses the fifth force 
on large scales; on small scales the suppression can be 
compensated by the decreases of m, which makes e~ mr 
larger, and the combined effect can be a weakened flat- 
tening of AP/P again. 

3. Changes in r (and similarly in £) make either to or 
/3 more sensitively dependent on local matter density, 
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FIG. 20. (Colour online) The ratio between the mass functions of the dilaton models and the ACDM paradigm at a — 1.0. 



the deviation from linear perturbation results and the 
screening effect get stronger, especially at late times 
when structures have developed. This explains why at 
late times AP/P can decrease with time when varying 
r and £. 

The above results imply that the shape of the nonlinear mat- 
ter power spectra can be different in dilaton and other modi- 
fied gravity (e.g. chameleon) models. This will be studied in 
more details in a forthcoming work. 



2. Mass functions 

This subsection contains the result of the mass functions 
from the dilaton simulations. The method to calculate the av- 
erages and standard deviations here is the same as that used in 
the symmetron case. 

Fig.l20lshows the results at a = 1.0, where we can see that 

1 . The dilatonic fifth force enhances the formation of dark 
matter structures. The effect is strongest for medium- 
sized halos and is weaker for very large and very small 
halos. As in the symmetron case, this is because for very 



large halos the screening effect weakens this enhance- 
ment, and many of the small halos have accreted more 
matter or merged with other halos to form larger halos. 

2. As discussed in § 1111 B 21 decreasing A 2 makes the fifth 
force less screened, and as a result more large halos are 
formed and fewer small halos survive the mergers . 

3. Increasing /3q makes the fifth force stronger and pro- 
duces more halos of all mass ranges probed by our sim- 
ulations. The dependence on /3q is quite sensitive, for 
example, for /3q = 1 the deviation from ACDM can be 
up to 50%, while for /3 = 0.25 this is less than 5%. 

4. As in the case of the matter power spectrum, the mass 
function becomes larger as r increases, and the depen- 
dence on r is quite weak, especially when r < 1 (mod- 
els C2, C3 and C4). As mentioned above, this is be- 
cause increasing r simultaneously increases the cou- 
pling strength and decreases the range of the fifth force, 
and the two effects cancel to some extent. 

5. The ^-dependence of the mass function shows a similar 
behaviour to that of the matter power spectrum. For ha- 
los more massive than ~5x 1O 13 /i -1 M , we find that 
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decreasing £ results in more halos being produced, sim- 
ilarly to the matter clustering power at k ~ 1/iMpc -1 . 
For smaller halos, model D2 predicts fewest while D3, 
D4 gradually catch up Dl, which is similar to the matter 
power at k > 3 — 4AMpc _1 . Overall, the ^-dependence 
is quite weak, similar to the r-dependence. 

As in the case of the matter power spectra, we are inter- 
ested in the shapes of the mass functions. As discussed above, 
changing r (or £) makes either the scalar field mass or the cou- 
pling strength more sensitively depend on local matter density, 
and in both cases the screening gets stronger (especially for 
large halos), consistent with what is seen in the matter power 
spectrum. A change in A2 strengthens or weakens the screen- 
ing effect but does not change the coupling strength for un- 
screened particles, and as a result the mass function behaves 
as in f(R) gravity models [28]. Finally, a change in flo mainly 
affects the coupling strength for unscreened particles, but not 
so much the degree of screening, which is why An/n flattens 
for large halo masses. 

To see how the dilaton effect on the mass function changes 
with time, we also show in Fig.|2T]the ratio between the mass 
functions at a = 0.5. As discussed in the previous subsection, 
at this time the linear perturbation theory is a better approxi- 



mation to the full theory. This implies that the screening of the 
fifth force has not yet been very significant, as is confirmed by 
this figure, which shows a weaker suppression of the dilaton- 
to-ACDM ratio at the high mass end. As in Fig.|20l the mass 
function results at a = 0.5 show a good match with the be- 
haviour of the matter power. Note also that the effect of vary- 
ing r and £ is larger at early times, which also agrees with the 
behaviour of matter power spectra. 

The above results indicate that the period between a = 0.5 
and a = 1.0 is an important era for the dilaton model, during 
which the structure formation is significantly affected by the 
nonlinearity of the model. In particular, we see that the shape 
of AP/P and An/n experiences qualitative changes during 
this period. 



VII. DISCUSSIONS, SUMMARY AND CONCLUSIONS 

A. Symmetron and dilaton screening 

Modified gravity models vary according to their screening 
mechanisms by which the fifth force is suppressed in local en- 
vironments. The Vainshtein mechanism works in theories of 



27 



the Galileon type where a scalar field with non-canonical ki- 
netic terms couples to matter in a reduced fashion in dense 
environments. Chameleons have an environment-dependent 
mass that becomes large enough to Yukawa suppress the fifth 
force in dense regions. Finally, the symmetron and the dila- 
ton share a similar mechanism whereby the coupling of the 
scalar field to matter is field-dependent and can vanish in the 
presence of dense matter. What distinguishes these two types 
of models is their scalar potentials: a Mexican-hat for sym- 
metrons and a monotonic function for dilatons. The coupling 
function for both types of models is a quadratic functioio 

Following the idea of Il6[[l7ll . the generalised dilaton and 
symmetron models studied here are completely specified by 
two temporal functions m(a) and /3(a). These give the most 
general models with a quadratic coupling to matter and scalar 
field mass that is a power-law function of a in the background 
cosmology for the generalised dilatons. For the generalised 
symmetron models, the scalar field mass vanishes for a < a* 
and increases to its present cosmological value from then. In 
both models, the screening of the fifth force is achieved in high 
density regions where the scalar field is trapped near the min- 
imum of A(tp). Yet the temporal dependences of the coupling 
to matter are drastically different: for generalised symmetrons 
it varies smoothly from a vanishing value for a < a» to its 
present value whereas the generalised dilatons it grows expo- 
nentially fast in the recent past of the Universe to reach its 
present value. 

As discussed in lfl7ll . the background expansion rate of such 
models is practically indistinguishable from that of the stan- 
dard ACDM paradigm, so that the cosmological effects of the 
fifth force could only be seen in the large-scale structures. In 
this work, we have performed large-scale A-body simulations 
for the generalised dilatons and symmetrons, investigating in 
detail the effects of varying the dilaton and symmetron pa- 
rameters on the nonlinear structures of the Universe. Some of 
these parameters are associated with the coupling to matter (3q 
(/3* for the symmetron case), and £ which specifies the range 
of the fifth force on the cosmological background. A few ex- 
tra parameters are used in the parameterisation to define the 
shapes of the potential and coupling function as functions of 
the scalar field. For the dilatons, these parameters are A2,r 
and for the symmetrons they are a», N and M. 

Let us first discuss the common features of these models: 

1 . The coupling to matter /3q (or /3*) determines the overall 
strength of the fifth forces, and increasing them leads to 
more structures. 

2. Decreasing £ leads to a shorter range for the fifth force 
and therefore a smaller enhancement of matter cluster- 
in£B 

In the end, the effects on structure formation are mainly de- 
termined by how fast the fifth force evolves and how efficient 



Of course, other types of coupling functions can be used, as we have done 
in the generalised symmetron model. 

In the dilaton case, changing § also affects the coupling strength, making 
the dependence on £ more complicated. 



it is screened in dense regions. An intuitive way to see this is 
to look at the expressions of /3(a) in these two models, as our 
discussion on tomography shows that this could be translated 
into (3(p m ), therefore giving us a sense about the screening, at 
least qualitatively. From Eqs. ( T53ll62b we can see that 

1. In symmetron models, the coupling vanishes at a < a* 
(or equivalently for p > p*) and after that it grows as 
a power-law function. Varying from to /3* between 
a = a* and today, j3 depends quite sensitively on a or 
p m in the regime with p m < p* ; however, the symmetry 
of V C H can be quickly restored for p m > p* resulting in 
a strong suppression of the fifth force. In other words, 
there is a clear cutoff density beyond which the screen- 
ing is very effective, and this cutoff is close to p*, which 
is fairly low. 

2. In dilaton models, the coupling grows exponentially 
with time and with decreasing density. As can be seen 
in Eq. d62b . /3 decreases and becomes vanishingly small 
if one goes back in time or goes to high-density regions, 
much more quickly than it does in the symmetron mod- 
els [c.f. Eq. (l53bl. This implies that the dilaton screen- 
ing can become effective for lower densities than the 
symmetron mechanism. 

It appears that the dilaton screening mechanism is more ef- 
ficient than the symmetron mechanism. However, local tests 
of gravity are carried out in very dense regions, where the 
fifth force can be strongly suppressed in both models. With- 
out specifying the exact parameter values for a given model, 
being it dilaton or symmetron, it is hard to say which one can 
satisfy local constraints more easiljOB 

B. Summary of numerical results 

Let us now summarise the results for each model. 



1. Generalised symmetron models 

The symmetron models we have simulated are close to what 
is allowed by local gravity experiments. Those constraints are 
mainly on the combination of the parameters a* and £ with the 
coupling strength /3* being an (almost) unconstrained param- 
eter. This parameter, which controls the magnitude of the fifth 
force compared with gravity, can in principle be constrained 
by its effect on the cosmic structure formation. 

Our simulations show that for a fiducial value of /3* = 1.0 
the symmetron models predict an enhancement of the nonlin- 
ear power spectrum with respect to ACDM of up to 40% for 
k ~ 1 /iMpc -1 and up to 50% at k ~ 10 /iMpc -1 . Likewise 
we find an enhancement of up to 50% in the mass function for 
halo masses in the range of 10 12 — 1O 14 /i _1 M0. 



It is clear that by varying the parameter values both models can be made 
either more or less screened. 
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We have shown how the fifth-force effect is changed by 
varying the other four model parameters: a*, iV, M and £. 

1. The parameter a» controls when the symmetry in 
V e ft(<p) is broken so the fifth force becomes non- 
vanishing. Decreasing a* gives it more time to influ- 
ence the matter clustering, as a result not only the mat- 
ter power spectra and mass functions deviate more from 
the ACDM results but also their shapes change qualita- 
tively (more discussion below). 

2. N is the parameter which controls the coupling strength 
via /3 oc <p N . Since \<p\ is very small, increasing N will 
suppress the magnitude of f3 (or the fifth force), and 
therefore causes less clustering of matter. 

3. M is the shape parameter of the symmetron field poten- 
tial, which determines how easy it is for if to roll away 
from ip — where f3 vanishes. Increasing M makes 
this easier, leading to a less-screened fifth force and thus 
more clustering and structures of matter. 

4. £ controls the scalar field mass and therefore the range 
of the fifth force in vacuum, A* = 2998£/i _1 Mpc. In- 
creasing £ makes the scalar field mass (range of the fifth 
force) proportionally larger (shorter), and thus leads to 
a stronger suppression of the fifth force and limits its 
range. 

As a rough guidance, increasing the symmetry-breaking 
scale factor a* from 0.33 to 0.50, decreasing A* from 
2.0/i _1 Mpc to 1.0/i _1 Mpc, increasing N from 2 to 4 or re- 
ducing AI from 6 to 4 are found to lower the enhancement of 
the power spectra and mass functions by ~ 10 — 20%. The 
parameters we adopt in the simulations are in the 'realistic' 
range and can be tested by future galaxy surveys. 

2. Generalised dilaton models 

We have also studied how structure formation in the gen- 
eralised dilaton models is affected by varying the four model 
parameters A%, /3q, r and £. 

1 . The effect of increasing is to make the total effec- 
tive dilaton potential V e s (if) steeper and so to keep the 
scalar field closer to <p*, where /3 and the fifth force van- 
ishes. The ACDM limit is retrieved by letting A2 —> 00. 
According to our simulations, reducing A2 to 5 x 10 4 
produces a ~ 20% enhancement in the nonlinear mat- 
ter power spectrum between z = 1 and z = 0, which is 
significantly smaller than the linear perturbation predic- 
tions, demonstrating the efficiency of the dilaton screen- 
ing mechanism. It also enhances the mass function by 
maximally ~ 25% in the same redshifts. These numbers 
assume that /?o = 0.5. 

2. The effects of increasing /?o are to strengthen the fifth 
force overall, and /3q = corresponds to the ACDM 
paradigm. The simulations show that even increasing 
/3o to 1.0 only causes 30 — 35% enhancement in the 



matter power for scales smaller than k ~ 1/iMpc -1 be- 
tween z = 1 and z = 0. This is at least 50% smaller 
than the linear perturbation result, again showing that 
the fifth force is efficiently screened in dense regions. 
In the mean time, the mass functions are increased by 
up to 50% with respect to the ACDM prediction. These 
numbers assume that A2 = 10 5 . 

3. Increasing r to 3/2 simultaneously increases the 
strength and decreases the range of the fifth force. The 
r-dependence of the matter clustering is rather weak as 
a result of the cancellation due to these two opposite 
effects. Assuming A 2 = 10 5 and /3 = 0.5, increas- 
ing r to 1.333 only enhances the matter power spec- 
tra by less than 10% at k — 1/iMpc -1 and 15% at 
k ~ 10/iMpc -1 , which is again significantly smaller 
than the predictions of linear perturbation theory. The 
mass function increases by up to 25% in this case. 

4. The effects of increasing £ are similar to those of de- 
creasing 7-, and as a result the dependence on £ is also 
fairly weak. 

Again, future galaxy surveys can place realistic constraints 
on the models studied here. 



3. Highlights and comparisons 

In both the generalised symmeton and dilaton models, as 
in f(R) gravity models l34Tl .we find that at late times the lin- 
ear perturbation theory fails to be a good approximation even 
for quite large scales (k ~ 0.05/iMpc -1 ). However, at earlier 
times it gives better agreement with the full simulations. This 
indicates that the environmental suppression of the fifth force 
becomes more important at late times when cosmic structures 
(very dense matter clumps) have already formed. This high- 
lights the importance of numerical simulations in the study of 
(screened) modified gravity models. 

The deviations of matter power spectra and mass functions 
from ACDM in the symmetron and dilaton models are not di- 
rectly comparable, because they depend on the exact param- 
eter values used in each model. However, we can see that the 
shapes of AP/P and An/n can be very different in the two 
models, which is probably a consequence of the different be- 
haviour of the respective fifth forces. 

At early times, AP/ P increases with k in both models (see 
e.g., Figs.[T4"land[T9ll, similarly to what we see in f(R) gravity 
models 028113411 . Differences appear at late time when the fifth 
force has been in effect for long enough: 

1. For f(R) gravity models we see that AP/P develops a 
peak at k ~ 0(l)/iMpc _1 , and on even smaller scales 
it decreases with k. The peak comes from the enhanced 
matter clustering due to the fifth force acting between 
clusters, and the turnover on small scales is because 
(compared with ACDM result) on these scales the short- 
range fifth force still accelerates particles and prevents 
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them from further clusteringPl 

2. In the symmetron case, we also see the peak of AP/P 
at k ~ 0(l)ftMpc _1 , and on even smaller scales it goes 
up again. This seems to imply that the particle veloc- 
ity inside halos stops being enhanced after the screen- 
ing effect has kicked in (recall that the 'cutoff den- 
sity for screening is quite low here and that 'screening' 
here means a suppression of the amplitude, rather than 
range, of the fifth force), as a result of which the shape 
of AP/P on small scales is preserved since early times. 

3. In the dilaton models, no obvious peak of AP/P can 
be seen: the power spectrum seems to have flattened 
on scales smaller than k ~ 1/iMpc -1 . Such a flatten- 
ing in AP/P is expected in the linear perturbation re- 
sults for both the symmetron and dilaton models, as 
in the linear regime the time at which the fifth force 
becomes non-negligible is scale-independent below the 
scale TOq *. For symmetrons the flattening is destroyed 
by the screening effect, while for dilatons it is not. 
As mentioned in § I VII Al dilaton screening can apply 
to lower matter densities: this indicates that the inter- 
cluster fifth force can be strongly suppressed as well, 
and thus the peak has not yet developed (notice that in 
some cases, such as B4, there is a small bump). Again, 
a more definite conclusion could only be drawn after a 
more detailed study of the density and velocity fields in 
the simulations, which is beyond the scope of this paper. 

The shape of An/n at late times is similar in symmetron, 
dilaton and f(R) gravity models, and the most important fea- 
ture is that it goes down in the high-mass end, demonstrating 
efficient screening of the fifth force in these large structures. 
At early times, however, An/n for the dilaton models show 
very weak mass dependence, which is close to the linear the- 
ory prediction, namely the fifth force is scale-independent. 

In the symmetron models, varying the parameters a^,N,M 
and £ changes the shape of AP/P (and of An/n) in similar 
ways, which results in a degeneracy in these parameters. This 
is because all these parameters control the degree of screening 
of the fifth force. 

This is not the case for the dilaton models, in which only a 
variation of Ai changes the screening mono tonic ally. Varying 
fto changes the overall strength of the fifth force more than its 
screening, while varying r or £ changes the screening in more 
complicated ways. As a result there is no degeneracy in these 
parameters, except between r and £ (see Fig.fTTIi. 



C. Conclusions and outlook 



1. to show the power of the modified gravity parameterisa- 
tion proposed in lfl6l [l7ll in systematic studies of struc- 
ture formation, 

2. to acquire a sense about the qualitative behaviour of the 
generalised symmetron and dilaton models, and the ef- 
fects of varying individual parameters, and 

3. to make a preliminary exploration of the 4-dimensional 
parameter spaces in these models and find models 
which are testable by the near-future observations. 

For all the test models in this paper, we find deviations from 
ACDM with similar magnitudes as those found in the f(R) 
gravity model J28ll34ll . which means that many of the cosmo- 
logical tests of f(R) gravity |29[ [3lU33ll could in principle be 
carried out here as well. 

On the other hand, the predictions of the cosmological ob- 
servables can be different from those in other modified grav- 
ity models with screening mechanisms, such as the chameleon 
models. For example, the shape of the matter power spectrum 
can be different in the symmetron, dilaton and f(R) gravity 
models, which implies that the respective screening mecha- 
nisms indeed work quite differently. It would be interesting 
to understand better the origin of such differences and see if 
they can be used to distinguish between the different modified 
gravity models in cosmology. These studies are under way. 
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In short, the aim of this paper is threefold: 



Contrary to intuitive understandings, this is not because 'the fifth force is 
suppressed on small scales'. The chameleon effect only reduces the range 
of the fifth force, but not its amplitude within that range. 
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